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We consider the Einstein-Maxwell system as a Cauchy initial value problem taking the electric 
and magnetic fields as independent variables. Maxwell's equations in curved spacetimes are derived 
in detail using a 3-1-1 formalism and their hyperbolic properties are analyzed, showing that the 
resulting system is symmetric hyperbolic. We also focus on the problem of finding initial data for 
multiple charged black holes assuming time-symmetric initial data and using a puncture-like method 
to solve the Hamiltonian and the Gauss constraints. We study the behavior of the resulting initial 
data families, and show that previous results in this direction can be obtained as particular cases of 
our approach. 
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I. INTRODUCTION 

One of the most important results in numerical relativ- 
ity in recent years has been the successful simulation of 
the coalescence of two spiraling (spin or spin-less) black 
holes (see Refs. [H, for an overview). This simulations 
are important since one expects black hole collisions to be 
among the most powerful sources of gravitational radia- 
tion. Gravitational radiation from this type of sources 
will presumably be measured by the next generation 
of interferometric gravitational observatories within the 
next decade or so It is therefore important to have 
simulations of different kind of astrophysical scenarios in 
order to compare with the observational results in or- 
der to reach a deeper understanding of such sources of 
gravitational radiation. In addition to the simulations of 
coalescing black holes, different authors have started to 
analyze collisions of extended objects (i.e. objects con- 
structed with a non-zero energy momentum tensor) like 
neutron stars 0,J^,JE 0, @| or even more exotic objects 
like boson stars [3, [ifll- Since such objects are less com- 
pact than black holes, and since their individual masses 
are limited, one expects that the amount of gravitational 
radiation emitted by the collision of these objects will 
be weaker relative to the two-black hole problem. How- 
ever, from the numerical point of view, such scenarios 
are by far more challenging since hydrodynamics, micro- 
physics or field theory are also involved. Moreover, their 
gravitational-wave signals can be richer in the sense that 
they can carry information about the internal composi- 
tion {e.g. equation of state). 

Another interesting numerical problem that one can 
conceive and that might have an observational counter- 
part is the collision of two charged black holes (TCBH). 
This problem is simpler than the case of two neutron stars 



•Electronic address: 'malcubi@nucleares.unam .mxl 
^Electronic address: jcdegollado@nucleares.unam.mx] 
■f Electronic address: ,imarcelo@nucleares.unam.mxl 



but perhaps more interesting than the two uncharged 
black hole collision. In fact, it is conceivable that if black 
holes in binary systems were formed by the gravitational 
core collapse of neutron stars or supernovae then they 
could have a small amount of charge. Actually, even if 
a black hole is originally uncharged but immersed in a 
uniform magnetic field (which in turn can be produced 
by an accreting plasma surrounding the back hole), it 
can be charged up to some extent Very likely their 
charge would be small compared with their mass (in suit- 
able units), but perhaps large enough to leave an im- 
print in the wave forms of gravitational radiation dur- 
ing a collision. In fact, numerical simulations of electro- 
magnetic fields immersed in the background spacetime 
corresponding to the collision of two uncharged black- 
holes already show that the dynamics of the background 
spacetime induce the emission of electromagnetic radia- 
tion that is correlated in a very particular way with the 
gravitational- wave signals [l^ . Such electromagnetic ra- 
diation (if detected) might provide information about the 
pre-merger stage as a precursor to the coalescence of the 
black holes. One can expect that taking into account the 
back-reaction of the electromagnetic field in the space- 
time itself can have even more interesting features. 

But even from the theoretical point of view, it seems 
important to analyze the interplay between the gravi- 
tational and electromagnetic forces in a TCBH collision. 
Indeed, the analysis of the interplay between electromag- 
netism and gravity within a BH spacetime has a long 
history in general relativity [i^. The first analytic so- 
lution to Einstein's equations involving an electro mag - 
netic field was given by Reissner and Nordstrom [l^, [l^ . 
This solution was interpreted as a spacetime containing 
a static and spherically symmetric charged black hole. 
Much more later, Papapetrou and Majumdar [l^ 
found a static solution involving multiple black holes hav- 
ing the maximal charge-to- mass ratios (see Sec. IIIIB[) . 
Perjes jlTj, and Israel and Wilson [l^ generalized this 
solution to the stationary case. Several uniqueness the- 
orems within the Einstein-Maxwell theory have been es- 
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tablished during the past which show how to character- 
ize certain kind of stationary black hole solutions (see 
Ref. [l^ for a review). 

More recently, Ruffini and collaborators ^2d\ have ana- 
lyzed the case of a static charge in a Reissner-Nordstrom 
spacetime by using a perturbative approach. In the 
framework of the so-called membrane paradigm [2]| , and 
using a 3-1-1 formalism, Thorne and colleagues have dealt 
with different problems involving electromagnetic fields 
in strong gravitational background fields (namely those 
generated by stationary black holes). Such a "mem- 
brane" viewpoint allows one to approximate several re- 
sults concerning electromagnetic fields around BH's, no- 
tably at the horizons. The membrane viewpoint assigns 
to the horizon thermodynamic, mechanical and electric 
properties. In this direction one can mention the pio- 
neering work of Damour [l^l and Znajek [2^, who ana- 
lyzed the boundary conditions of electromagnetic fields 
at the horizons of a BH. Such boundary conditions can 
be thought of as arising from the physical properties of 
a fictitious membrane residing at the horizon or near the 
horizon (see [2l| for the introduction of the notion of 
"stretched horizon" which allows to approximate these 
boundary conditions near the true horizon) . Since in this 
viewpoint one neglects the back-reaction of the matter 
fields into the spacetime, this approximation will break 
down in situations where the self-gravity of the matter is 
important. It is in this regime where numerical relativity 
becomes crucial. 

In the context of multiple charged black holes, the work 
by Bowen [23| was one of the first to address the ini- 
tial data problem. He considered the case of zero initial 
magnetic field, but without imposing a moment of time- 
symmetry and without resorting to electromagnetic po- 
tentials. We will often make reference to Bowen's work 
in this paper. Previous to Bowen's work, Lindquist psj 
generated initial data for many stationary charged parti- 
cles by imposing time-symmetry and using a method of 
images and electromagnetic potentials. The present pa- 
per is similar in spirit to Bowen's work, except that we 
shall mainly be interested in a TCBH initial data which 
is computed using a method analogous to the puncture 
method p6j which is not inversion symmetric (see Sec. 

imi. 

In order to tackle the problem of a TCBH collision 
there are many challenges, both numerical and analyti- 
cal. First, from the numerical point of view there are two 
main considerations: one is the initial data and another 
one is the evolution of the Einstein-Maxwell system. For 
the former one requires suitable initial data compatible 
with the constraint equations. Whereas for the latter one 
needs a numerical code to solve the Einstein-Maxwell 
evolution system together with all the numerical tools 
necessary to analyze the location of horizons, the amount 
of emission of gravitational and electromagnetic radia- 
tion, etc. 

As regards the analytical challenges, there are various 
issues with different levels of complexity. Since the vast 



majority of numerical relativists are concerned with the 
3-1-1 formulation of Einstein equations and their corre- 
sponding numerical solution, the first step toward our 
goal consists in obtaining a well defined 3-1-1 decompo- 
sition of Maxwell's equations in a curved spacetime. It 
turns out that this problem has been embraced in the 
past in at least two works: first by Ellis in [131, and 
later by Thorne & Macdonald in [2^. In Section 
we present a completely independent derivation of the 
3-1-1 Maxwell equations which we then compare with the 
one reported by Thorne & Macdonald. Our derivation 
is based primary in the 3-1-1 formulation considered by 
York [29|] (see Refs. ^ lU for a thorough review of 
this formulation). As we will show, the 3-1-1 Maxwell 
equations in curved spacetime are in fact very similar 
to the usual Maxwell equations in flat spacetime, except 
for the fact that some extra terms due to the curvature 
arise. We obtain two "scalar" constraint equations for 
the electric and magnetic fields, and two "vector" evolu- 
tion equations for those fields. In this case, the electric 
and magnetic fields are referred to the so-called Eulerian 
observers whose four-velocity is orthogonal to the space- 
like hypersurfaces that define the foliation of the four 
dimensional spacetime (see Sec. [Tl]). The advantage of 
this 3-1-1 formulation of Maxwell equations, as opposed 
to the one based on electromagnetic potentials, is that 
their fundamental variables are gauge invariant ab-initio. 
Therefore, one needs to focus only in the gravitational 
gauge issue. Moreover, the set of equations turns out to 
be manifestly hyperbolic (symmetric hyperbolic). The 
only difference with respect to the flat case is that the 
eigenvalues and eigenvectors associated with the princi- 
pal part of the equations include contributions due to the 
(densitized) lapse and shift. We analyze these aspects in 
Sec. US 

As we mentioned above, another challenge which can 
be both numerical and analytical concerns the initial 
data. Let us recall that in the case of the two uncharged 
BH collision one can consider a rich family of interesting 
initial conditions. The simpler ones consist on impos- 
ing a moment of time symmetry in which the two BH 
are initially at rest and without angular momentum and 
spin. With such a condition the momentum constraints 
are satisfied trivially. In this case, the initial data will 
generate a head-on collision. Moreover, one can assume 
that the 3-metric is conformally flat. Both conditions 
in the Hamiltonian constraint lead to an elliptic equa- 
tion for the conformal factor (see Sec. IIII A|) . A unique 
solution of that equation is obtained when imposing suit- 
able boundary conditions which in turn are related to the 
topology of the initial hypersurface Sq (see Refs. [lO, HI] 
for a pedagogical review on initial data). One possible 
topology for Sq is R"^ minus two balls (the boundary of 
which represent the horizons of the two BH). This kind 
of initial data has been used in the so-called excised ap- 
proach where one ignores what is inside those balls while 
evolving the rest of the spacetime. A simpler topology 
for Eq is to consider M'^ minus two "points" . This is the 
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so-called "puncture data" approach. The punctures rep- 
resent in fact two different asymptotically flat regions on 
a two-throat shaped spacelike Sq. Since in vacuum the 
elliptic equation for the conformal factor is indeed linear 
one can in fact construct puncture data that represent 
an arbitrary number of black holes that are initially at 
rest. This data has a simple analytical expression whose 
form resemble the electric potential generated by a series 
of point charges. 

A much more realistic initial data was the one used 
in the simulations of the coalescence of two binary BH 
that we alluded above. In this case, there is a remark- 
able analytical solution (the Bowen-York solution) of the 
momentum constraints which represents two BH with ar- 
bitrary linear momentum and spin [s^. This solution is 
constructed using a conformal and transverse-traceless 
decomposition [32]. The most difficult part of this ap- 
proach is to solve a highly non-linear elliptic equation for 
the conformal factor. Again, one can adopt the puncture 
or the excised approach. In either case one needs to solve 
the elliptic equation numerically. The puncture approach 
has been very popular recently and is the one that has 
led to many successful evolutions of different BH binary 
configurations. The outcome of such evolutions has led 
to the prediction of a large variety of gravitational wave 
forms that in the near future will be confronted with the 
observational data Q . 

Returning now to the case of the two-charged BH prob- 
lem, one can take advantage of the experience with the 
uncharged case in order to construct interesting initial 
data. Unlike the uncharged case, the new difficulty is 
that the Einstein constraint equations will have contri- 
butions due to the electromagnetic fields. In order to sim- 
plify the problem the first attempt consists in assuming 
again a moment of time symmetry where the magnetic 
field is zero. As mentioned before, this kind of initial data 
was considered in the past by Bowen [23] ■ In this case 
the Poynting vector (which is the source of the momen- 
tum constraints) becomes identically zero. We could then 
use the Bowen-York initial data for this problem. How- 
ever, one still needs to solve the Hamiltonian constraint 
plus the Gauss constraint for the electric field, the for- 
mer containing now the contribution due to the electro- 
static energy-density associated with the initial electric 
field generated by the two-charged BH. One can further 
simplify the problem if one assumes that the two BH are 
initially at rest (zero spin and zero linear and angular 
momentum). This initial data would then represent the 
head-on collision of two charged BH. However, unlike the 
vacuum case, the elliptic equation for the conformal fac- 
tor is now highly non-linear. Remarkably, we have found 
a way to solve analytically both constraints (the Hamil- 
tonian and Gauss constraints) using a puncture approach 
(see Sec. IIIIB]) . This solution represents a superposition 
of multiple charged black holes all of which having the 
same charge-to-mass ratio. When this last condition is 
dropped, finding an analytical solution seems difficult. 
However, it is not difficult to find numerical solutions 



(see Sec lmC]) . 

The paper is organized as follows: Section HIl presents 
our derivation of the 3-1-1 Maxwell equations in a curved 
spacetime, as well as the analysis of their hyperbolic 
properties (this section is complemented by an Ap- 
pendix). The Einstein-Maxwell system is summarized 
and a brief discussion on the electromagnetic potentials 
is also included. In Section IIIII we analyze the initial 
data for multiple black holes using a conformal approach, 
where both analytical and numerical results are obtained. 
Finally Section IIVI contains several comments and re- 
marks for the future. 

II. THE MAXWELL EQUATIONS IN 3+1 FORM 

In the following we assume that the reader is famil- 
iar with the 3-1-1 formalism of general relativity (see 
Refs. d^, [13, [Slj, for a thorough review, and [sl] for 
the conventions adopted here). 

Let us first remember that the covariant Maxwell equa- 
tions read 

Va^^"' = -Anj" , (2.1) 

where 

Fab = daAb-dbAa, (2.3) 
p*ab _^^abcdp^^^ ^2.4) 

are the Faraday's tensor and its dual. Here we take the 
convention that e'^^^'^ = — 1/-^^ and 60123 = +^z^--3, 
with the signature of gab taken as (— , -|- , -I- , -I-) [43] . 

In order to obtain the 3+1 decomposition of Maxwell 
equations from their covariant form, on has to proceed 
in a way very similar to the derivation of the so called 
Arnowitt-Deser-Misner (ADM) equations of General Rel- 
ativity. Let us briefly review the ingredients of this de- 
composition. One considers a spacetime (Af , gab) (as- 
sumed to be globally hyperbolic) which is foliated by a 
family of spacelike hypersurfaces St {t e M) parametrized 
by a global time function t {i.e. M has topology 
M = St X R). The foliation is achieved in the following 
way: On a Cauchy surface Sf , one is given an initial data 
set that satisfies some constraint equations (see Sec. Ill Cp . 
The full spacetime is "reconstructed" by evolving these 
initial data using a suitable set of evolution equations 
(which includes the gauge) . Such a set of constraints and 
evolution equations are known as the ADM equations of 
General Relativity. 

In order to find a similar set of equations for the elec- 
tromagnetic case, an analogous algebraic and geometric 
decomposition of the field equations (|2.ip and (|2.2p has 
to be performed. The general procedure for obtaining 
a 3-1-1 splitting of a system of covariant field equations 
(also called orthogonal decomposition) consists on pro- 
jecting the different tensor fields in the directions par- 
allel and orthogonal to the timelike unit vector field n'^ 



4 



{rf'na = — 1) which is normal to St. The projection onto 
Tit is performed by first defining the projector operator. 



(2.5) 



This tensor field has the property of being idempotent: 

A tensor field ^r^^"' -"';^^^ is said to be tangen- 
tial to Sf if, when contracted with n'^ it gives zero, or 
equivalently if contracted with h\ it remains unchanged. 
For brevity, such tensors will be termed 3— tensors. Any 
tensor field can be decomposed orthogonally by using /i"^ 
and In particular, a 4- vector w°' is decomposed as 
follows 



(2.6) 



where w° :— h\w and w±_ := —Ucw'^. Moreover, the 
3+1 splitting of the metric reads 



-{N^ - N'Ni)dt^ -2Nidtdx' + h,j dx'dx^ , (2.7) 



where the lapse function iV > is defined as to normalize 
the (future pointing) dual vector field ria = —NVat- The 
shift vector is given by iV := — where t° = [d/dt)"" 
is a vector field that represents the "flow" of the time 
lines and which satisfies f^V at = 1. This means that 
is orthogonally decomposed as t°- = —N°- + Nn°-, with 
A'' = —Uat"". Here hij is the 3- metric (or induced metric) 
of the manifold St. In order to avoid confusion a note 
on notation is important at t his p oint: many numerical 
relativity references (see e.g. [23, [13]) use a and /?° to 
denote the lapse and shift instead of A^ and A°, with 
a — N and /J" = — A^" (do note the change in sign in the 
shift!). 

Another important object is the extrinsic curvature of 
the embeddings St which is defined as[50t 



(2.8) 



where Dahbc = 0. Finally, one should mention that the 
indices of 3-tensors can be raised and lowered with ft,"'' 
and hab, and also their contravariant time components 
are identically zero. Moreover the projector /i*^^ applied 
to any 3-tensor field acts as a ^'j,. 

The Maxwell equations written in the form of an 
initial-data (Cauchy) problem can be obtained by pro- 
jecting Eqs. (12.11) and (12. 2p orthogonally and tangentially 
to the space-like hypersurface St using and the projec- 
tor operator (12. 5p . respectively. Note that with respect 
to a local coordinate basis adapted to the foliation we 
have Ha = (-A^, 0,0,0) and n" = {l/N,NyN). In par- 
ticular, in a flat spacetime one has Ua — (—1, 0, 0, 0) and 
72" = (1,0,0,0). 

We will need to deflne suitable quantities associated 
with the electromagnetic (EM) fleld before proceeding 
with the 3-f 1 decomposition of the Maxwell equations. 
In order to do this, consider first the 3-1-1 decomposition 
of an arbitrary 2-0 tensor, say H""^, which is given as 
follows 



H 



_ {3)jjab {3} jj±b 



with 



{3)jjab 
(3)^_Lb 
{3)fja± 



H 



-_L_L 



jcd 



HaribH 



ab 



(2.12) 



(2.13) 
(2.14) 
(2.15) 
(2.16) 



For the particular case where H"^^ — F°\ the antisym- 
metry of F"'' impfies that F^^ = 0, so [Hj 



(2.17) 



where /!„ stands for the Lie derivative along n". From 
the above definition one can obtain the following identity 



K. 



ah 



(2.9) 



which shows that Kab is in fact a 3-tensor field. Further- 
more, its trace is given by 



We now define the electric and magnetic fields as mea- 
sured by the Eulerian observers with four velocity in 
the following way [s^ 



-ribF 



ba 



(3) p±a 



B" := -UbE*"" = 



ba ^ (3) 



E 



.±a 



(2.18) 
(2.19) 



K 



(2.10) 



As is well known, the set (Tt, hab, Kab) provides the 
initial data for the gravitational field. This data in fact 
cannot be chosen freely, and has to satisfy the Einstein 
constraint equations (see Sec. Ill Cl . 

At this point it is useful to introduce a covariant deriva- 
tive operator compatible with hab- Given a 3-tensor field 
'T"^""-\b,...,,,one defines 



h%...h'^,^X:...h,fh{Wf'T^ 



did-2 



The last equalities in Eqs. (P7T5)l and ((^4^ result by 
noticing that the above definitions and the antisymmetry 
of E'''' and E*'''' both imply that E'' and B° are in fact 
3- vector fields since riaE"" = = UaB"", or equivalently 
E" = h^-^E"^, B° = h'^^B'^. Otherwise such equalities 
can be corroborated simply by inserting Eq. (|2.17[) and 
its corresponding dual in the definitions above. In this 
way, one finds that 



E' 



(3)^afc^^a£;b_£;a^fc 



(2.20) 



(2.11) On the other hand, from Eqs. (EUD, and ^^ITM we 
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have 



Tja bacd tt' 



2 
1 

= ^me 



J)acd 



bacd{3)p^^^ (2.21) 



where in the last equahty we used the fact that the con- 
traction of the totaUy antisymmetric Levi-Civita symbol 
gbacd f^i^Q symmetric tensors rtbTT.^ and ni,nc vanishes 
identically. Moreover 



where 



(3)g±a/s 



(2.22) 



(2.23) 



The last identity arises when using Eq. (|2.5|) in the above 
definition, plus the fact that all the contractions between 



^becd 

clliu iliuic 

Now, since e°^^^ 



and more than one factor Ua vanish identically [52 



= -l/(^^v^), where 
h:=det{hij) and (^h-^^fa = ' - N h'' „hf M^e'^"'"^ . it is 
clear that (3)e-Li23 ^ 
simply identify 



c'- d^ 



0123 



1/Vh- Therefore we 



(3)ga/ff _ (3)g-La/g ^ (2.24) 

with the 3-Levi-Civita symbol defined in such a way that 
(3)gi23 _ \l\fh [s^l- It is also convenient to introduce a 
"flat" Levi-Civita symbol as 



abc 



^abc 



be 



1 



(3) 



^abc 1 



(2.25) 
(2.26) 



where e^*"^ and eff,^, take values (0, ±1) as in flat space. 
We can also invert Eq. (|2.22p as follows 



{3)pab _ {3)^abcj^^ ^ 

SO that Eq. (l2?20)l now reads ^ 

pab ^ {3)^abc^^ _^ ^apb _ pa^b 



(2.27) 



(2.28) 



In the same way we can obtain the following 3-f 1 de- 
composition of F*°-^: 



p*ab ^ _ i3)^abcp^ ^a^b _ ^a^ 



(2.29) 



We note that, just as in case of flat space, the dual op- 
eration maps the electric and magnetic fields as follows: 

Ea ^ Ba, Ba ^ -Ea M- 



We are now in the position of performing the 3-1-1 split- 
ting of Maxwell equations (details can be found in the 
Appendix). The projection of Eqs. ()2.ip — ()2.2p onto ria, 
after the use of Eqs. (|2.28p and (|2.29p . leads to the ini- 
tial value constraints for the electric and magnetic fields 
respectively: 



DaE"- = 4^p, 
DaB'' = 0. 



(2.30) 
(2.31) 



where we remind the reader that Da is the deriva- 
tive operator compatible with hab [cf. Eq. (|2.1ip ]. and 
p := —Haj"' is the charge density as measured by the 
Eulerian observer. More specifically, the covariant 3- 
divergence is given by 

DaE'' ^^d,(^VhE'^ . (2.32) 

On the other hand, the projection of (|2.ip onto St 
provides the evolution equation for the electric field (see 
Appendix for details): 

h\£nE^ - 

{D X B)" - {By. af + KE" - An f , (2.33) 



with 



{D X BY 
(B X a)° 



(3)^.a 



(2.34) 
(2.35) 
(2.36) 



and where K :~ K\ is the trace of the extrinsic curvature 
given by Eq. ^T^, and a" := n^Van" = D%lnN) is the 
acceleration of the Eulerian observer (29l . [3]| . 

Taking the spatial components of Eq. (|2.33p one 
finds |53| 



dtE^ + Cf<iE' = {Dx NBY 

+ NKE' - Attn ^^'^f 



(2.37) 



where now £n is the Lie derivative along the shift, and 
where we used the fact that {D x NBy = N{D x Bf - 
N{B X a)\ 

In a similar way, the projection of (j2.2p onto Et pro- 
vides the evolution equation for the magnetic field 



dtB' + Ct<sB' ^-{Dx NEf + NKB' 



(2.38) 



A self-consistency check of Eq. (|2.37l) can be performed 
by noting that Eq. (|2.ip implies the charge conservation 
equation 







(2.39) 



This equation can be written in terms of the 3-1-1 lan- 
guage as follows (see Appendix): 



dtp + C^p = -Da{N (3)j-) + NpK 



(2.40) 
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Therefore, by replacing in the above equation the values 
of p and ('^^j'^ given by the Eqs. (|2.30p and (|2.37p respec- 
tively, one finds (after some algebra that involves the use 
of the commutator of covariant derivatives applied to a 
vector field as well as the Gauss-Codazzi equations) an 
identity. In a similar way, one can also check the self- 
consistency of Eqs. ([OT|) and 1^^ . 

At this point it is important to mention that a set of 
3-1-1 Maxwell equations analogous to Eqs. (|2.30p . (I2.3ip . 
(|2.37p and (|2.38p . as well as the 3-1-1 charge conservation 
Eq. (|2.40p ^were derived previously by Thorne & Mac- 
donald in [2^ u sing the same sign conventions but a dif- 
ferent notation [5a |. These authors in turn used the 3-1-1 
congruence formalism of Maxwell's equations derived by 
Ellis [13. 

It is important to note that in a flat spacetime 
Eqs. ((OO)) . ((OT|) . (p:T71) and reduce to the famil- 

iar form of Maxwell's equations. For the present case of a 
curved spacetime, one can in fact rewrite equation (j2.30p 
in integral form as 

/ E^VhcTada I -Uafy/hdx^dx^dx^. (2.41) 

On the left hand side (l.h.s), one has the flux of the elec- 
tric field lines across a closed two-surface lying on with 
normal aa G T^*. On the right hand side (r.h.s), one has 
the total charge measured by the Eulerian observers con- 
tained in the volume enclosed by the two-surface. The 
r.h.s is a consequence of Eq. (|2.39p . which implies that 
the total electric charge Q is conserved: 

Q ^ [ Nf\fh,dx^dx^dx^ 

= / -Uaj'^Vh dx^dx'^dx^ 



p\fh dx^dx^dx^ 



(2.42) 



St 



Note that in all the above expressions there appears the 
proper volume clement ^/hdx^dx'^dx^ on St as measured 
by the Eulerian observers. A similar analysis can be done 
on Eq. (|2.3ip . except that in this case there are no mag- 
netic charges (the magnetic field lines are always closed). 

As concerns the evolution equations (|2.37p and (|2.38p . 
some of the extra terms appearing there are due to curva- 
ture effects plus the fact that all observables are referred 
to the Eulerian observers. For instance, if one uses dur- 
ing the evolution the so-called maximal slicing condition 
(which is defined by imposing K = = dtK; this leads to 
an elliptic equation for the lapse N [cf. Eq. (|2.8ip ]). then 
the terms proportional to NK on the r.h.s of Eqs. (|2.37p 
and (|2.38p vanish identically. Otherwise, those terms are 
present and are associated with the time variation of the 
proper volume elements on S(. Now, quite independently 
of the choice of a particular time slicing, Thorne & Mac- 
donald have provided geometrical interpretations of the 
extra terms that couple gravity with electromagnetism 



in a non-trivial fashion. Such interpretations can be- 
come even clearer when writing the evolution equations 
in integral form [l^ . 



Hyperbolicity analysis of Maxwell's equations 
in curved spacetimes 



The system of evolution equations (|2.37p and (|2.38p 
for the electric and magnetic fields can clearly be written 
as [Ei 



diU 



(2.43) 



where u = (£^', B') are the fundamental variables and M* 
are the characteristic matrices along the directions x*. 

In order to gain some insight on the hyperbolic struc- 
ture of Maxwell's equations, let us first focus on the case 
of a flat spacetime background in Cartesian coordinates 
and in vacuum (i.e. in the absence of electric charges 
and currents). In such a case the matrix is a 6 x 6 
block antidiagonal matrix which can be written as the 
following direct sum of two 3x3 matrices: 



low J 



(2.44) 



where the symbol © means a direct sum by placing the 
blocks in the antidiagonal: 







low 



nip 



(2.45) 



The components of the upper and lower matrices 



and 



low 



are given respectively by 



ilm 



= e 



qiln 

gihn ih 

how — 



(1 < Z,TO < 3) , 
' (1 < ;,m < 3) 



(2.46) 
(2.47) 



The minus "— " in Eq. p.47p corresponds to the asymme- 
try in sign in the dynamic Maxwell equations (|2.37p and 
d^nS]). However, since Mjf^ = M™' this shows that 
is indeed symmetric [e.g. see Eq. (|2.52p ]. Therefore one 
concludes that the system of equations (|2.37p and (|2.38p 
are in fact symmetric hyperbolic and so they admit a well 
posed Cauchy problem. 

One can further analyze the different modes of the 
characteristic matrix M*. Let us take first the simple 
case where and _B* consists of plane waves moving in 
the 'x' direction so that 



(2.48) 
(2.49) 



Equations (|2.37p and (|2.38p then reduce to the following 
algebraic system (remember that we are considering flat 
spacetime in vacuum) 



ujE' + e'^^kj& = 
uI3' - e'^^kiE^ = 



(2.50) 
(2.51) 
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The matrix then takes the form 



/O 0\ 
c 
-c 

-c 
VOcOOOO/ 



(2.52) 



where c := k/u corresponds to the speed of hght. The 
constraints equations (|2.30p and (|2.3ip impose the fol- 
lowing conditions 



Ex = Q — Bx 



(2.53) 



Now, the eigenvalues of are A — (0,0, ±c,±c). 
The eigenvalue A = corresponds to the eigenvec- 
tors ei = (1,0,0,0,0,0) and e2 = (0,0,0,1,0,0), 
the eigenvalue \ — —c corresponds to the eigenvec- 
tors 63^(0,-1,0,0,0,1) and = (0,0,1,0,1,0), and 
the eigenvalue A — c corresponds to the eigenvectors 
eg = (0,1,0,0,0,1) and eg = (0,0,-1,0,1,0). 

The determinant of the eigenvector matrix is 



det (I 



-4 . 



(2.54) 



The set of eigenvectors therefore is complete, and the 
evolution system turns to be strongly hyperbolic. This of 
course is not a surprise since we already knew that is 
symmetric and therefore the system is in fact symmetric 
hyperbolic. 

The first two modes with speed zero (A = 0) are clearly 
unphysical since they are associated with modes that vio- 
late the constraints (|2.53[) . On the other hand, the modes 
with speed ±c do satisfy the constraints. These physi- 
cal modes correspond to 63, 64, es, eg, and are associated 
with the two polarizations states (each one with speed 
±c) of the electromagnetic waves. 

Now, for the case of waves propagating along an arbi- 
trary direction defined by the unit vector s, the principal 
symbol of the system is given by C := WVsi. Then we 
have 



/ 






S3 

V -S2 Si 







-S3 S2 

-Si 






-S3 

S2 







S3 

-Si 






-S2\ 
Sl 





/ 



(2.55) 



The eigenvalues are A = (0, 0,±1,±1) with correspond- 
ing eigenvectors (from now on we will take the speed of 
light to be equal to 1): 



• A = 0: 



ei = (si/s3,S2/s3,l,0,0,0) , (2.56) 

62 = (0,0,0,S1/S3,S2/S3,1) . (2.57) 



• A = -1 

6*3 = (s2,-(si -|-S3)/S1,S2S3/S1,-S3/S1,0, 1), (2.58) 
6*4 = (-S3,-S2S3/si,(Si +S2)/Si,-S2/Sl,l,0). (2.59) 



• A 



66 



(-S2, {sl + s2)/si, -S2S3/S1, -S3/S1, 0, 1). (2.60) 
(s3, S2S3/S1, -(si + S2)/si, -S2/S1, 1, 0). (2.61) 



and the determinant of the eigenvector matrix R is 



det 



'^l{44) 



(2.62) 



The set of eigenvectors is clearly complete, and the evo- 
lution system again turns out to be strongly hyperbolic. 
Since C is symmetric, the system is in fact symmetric 
hyperbolic. 

Again, the two modes with zero speed violate the con- 
straints since S ■ ei and S* • 62 do not vanish, where 
S :— (si, S2, S3, Sl, S2, S3) (the constraints correspond to 
^^u' ^ 0, that is, u''S^ = = E's, + B'si). On the 
other hand, one can see that S ■ Ci vanish identically 
for i = (3, ...,6). Therefore such modes satisfy the con- 
straints and propagate with the speed of light. They 
correspond to the two polarization states with speed ±1. 

We can now consider the full system of evolution equa- 
tions (|2.37p and (|2.38p with a prescribed shift. The prin- 
cipal symbol now reads, 



/ 


Ns 











AAs3 


-AAS2 \ 







Ns 





-AAS3 





JVsi 










Ns 


AAS2 


-AAsi 










-Ms; 


Afs2 


Ns 










AAs3 





-Msi 





Ns 





V 


-Afs2 


Nsi 











Ns J 



(2.63) 

where Ns :— N^Si, and with Af := N/Vh the densitized 
lapse. 

The eigenvalues of the principal symbol are now 
A = {Ns, Ns, Ns iA/"). These correspond to the charac- 
teristic speeds as measured by the Eulerian observers. In 
the flat spacetime limit they reduce to those previously 
obtained. Their corresponding eigenvectors are now 

• \ = Ns 

ei = (si/s3,S2/s3,l,0,0,0) , (2.64) 

62 = (0,0,0,S1/S3,S2/S3,1) . (2.65) 

• X = Ns-Af 

6*3 = (s2,-(Si +S3)/si,S2S3/si,-S3/si,0,l). (2.66) 
64 = (-S3,-S2S3/si,(Si -|-S2)/Si,-S2/Si,l,0). (2.67) 

• X = Ns+M 

6*5 = (-S2,(s?-HS3)/si,-S2S3/si,-S3/si,0, 1), (2.68) 
6*6 = (s3,S2S3/si,-(Si-|-S2)/si,-S2/si,l,0). (2.69) 
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Just as in the flat case, the first two modes which prop- 
agate with the coordinate speed Ng are unphysical since 
they violate the constraints. On the other hand, the re- 
maining modes satisfy the constraints and propagate at 
the speed of light ± Af. 

The characteristic speeds along a given direction can 
also be written as A' = {N\ N\ ± J\fs'), so that the 
projection along the s direction provides A. This shows 
that physical modes propagate along the light cones, cor- 
responding to the coordinate speeds ± Afs^ . 



B. The energy- momentum tensor of the 
electromagnetic field 



by 



The energy-momentum tensor of the EM field is given 



47r 



FacFf^' 



1 



:gab FcdF' 



cd 



(2.70) 



Using Eq. ([2?28l) we obtain 



FacF^ ' 



— {EaEb + BaBb) 

B^hat + E^nanb + 2E'B''^^h,d(anb} , (2.71) 



where E'^ = E'^Ea and = B'^Ba- 

From the above equation one can easily find 

FacF"" = -2 {E^ - S^) . (2.72) 

The energy-momentum tensor (|2.70p then becomes 



1 

47r 
1 



-{EaEb + BaBb) + -hab{E^ + B^) 



-naUbiE' + B') + 2E'B^ ^^^e.^ianb) 



, (2.73) 



where we have used the fact that gab = hab — naUb [cf. 
Eq. (|2.5p ]. The 3+1 decomposition of this tensor is [cf. 
Eq. (irT2)l ] 



Tab = ^naUb + UaJb + J a^b + Sab 

where now [cf. Eqs. (PT^ - ipH]) ] 
1 



b^—iE' + B') 
on 



(2.74) 

(2.75) 
(2.76) 



5*06 hahb'^Tcd 

= ^ [hahi.E'' + B^) - 2{EaEb + BaBb)] . (2.77) 

We identify £ with the energy-density of the EM field as 
measured by the Eulerian observers, Ja with the momen- 
tum density measured by those observers (the Poynting 
vector), and Sab with the stress tensor (60j . 



Since the trace of the energy-momentum ten- 
sor (I^TTO)) vanishes, then Eq. (P?7i)) leads to [cf. Eqs. 
(12751) and (j2T7)) ]. 

S = £, (2.78) 
where S — S'^a is the trace of the 3-tensor (|2.77p . 

C. The Einstein-Maxwell system 

We now consider the Einstein equations in 3+1 form, 
with the matter sources provided by the electromagnetic 
field contributions of Sec. Ill Bl 

The Hamiltonian and momentum constraints are, re- 
spectively [n, M, mi 



167r£ 



(2.79) 



- D'K = 8nJ' . (2.80) 
The dynamic Einstein equations read [2^ [13, l3lj | 

D,D,N 



dtKij 



KK,, 



2KaK. 



47riV [h,,{S^£)-2S,,] 



inNSi. 



(2.81) 



where we used Eq. (|2.78p to simplify the r.h.s of 
Eq. (|2.81|) . Taking now the trace in Eq. (|2.8ip . and using 
Eq. (|2.79p . one obtains the following evolution equation 
which can be very useful in many cases {e.g. see Sec. IIIII 
below) : 

dtK + N^diK + D'^N - NKijK^' 

= 47r7V (S" + £) = 87riV£ , (2.82) 

where we used Eq. (|2.78p in the last step. Here stands 
for the Laplacian operator compatible with the 3-metric 

hab- 

In all the four equations (|2.79p - (|2.82p . the r.h.s is given 
in terms of the energy-momentum contributions defined 
in Eqs. (|275l) - ([2J7p . 

It is well known that the evolution equations (|2.8ip . 
when written in first order form, are only weakly hyper- 
bolic (see Ref. [soj for a thorough review) , and so do not 
admit a well-posed Cauchy problem (in the Hadamard 
sense). However, by adding suitable multiples of the con- 
straints and using a conformal decomposition {e.g. the 
BSSN formulation [3^) one can write the evolution sys- 
tem in such a way that the evolution system admits a well 
posed Cauchy problem. Here we are not concerned with 
with that issue, but will will consider it when we study 
numerical evolutions of the Einstein-Maxwell system in 
a future paper. 

The Cauchy problem in this case can be summarized as 
follows: given the initial data (St, hab, Kab, Ea, Ba), sat- 
isfying the constraints ^1^, dMO]), (j^:^ and ((OT|) . 
one can then evolve forward in time (given a prescrip- 
tion for N and TV*) the fields hab, Kab, Ea, Ba using their 
evolution equations (|2.8p . (|2.8ip (or their correspond- 
ing strongly hyperbolic formulated equations), (j2.37p and 
(|2.38p respectively |,61i] . 
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D. Electromagnetic potentials 

Up until this point we have worked directly with the 
electric and magnetic fields and ignored the potentials. 
This is quite dehberate, as the use of the potentials can 
complicate matters, in particular due to the fact that one 
needs to choose a gauge. Also, the evolution equations 
for the electromagnetic field when written in terms of 
the potentials are second order in space and time, which 
brings extra complications coming from the fact that co- 
variant derivatives do not commute. Nevertheless, here 
we will very briefly describe, without going into any de- 
tails, the 3-1-1 form of the electromagnetic potentials and 
their relation with the electric and magnetic fields. 

Let us start by remembering that in terms of the po- 
tential 4-vector A° the Faraday tensor is given by equa- 
tion (|2.3p . which we rewrite here for concreteness: 

Fab = daA, - dbAa . (2.83) 

Starting from the potential 4-vector, one can now define 
a 3-1-1 "scalar potential" $ through 

$ := -UaA'^ , (2.84) 

together with a potential 3-vector (3)^4° defined as 

(3)^a ._ i^a^^b (2.85) 

From these definitions one can inmediately find that, in 
a coordinate system adapted to the 3-1-1 foliation: 

$ - iVA* - -1 {At + N^Aa) , (2.86) 

and 

Now, by projecting the expression for the Faraday ten- 
sor in terms of A'^ given above one can obtain, after some 
algebra, the following relation between the electric-field 
Ei and the 3-f 1 potentials: 

dt (^U, + £n ^^^A, - -NE, - A (iV$) . (2.88) 

Notice that this equation can in fact be interpreted as an 
evolution equation for the potential 3-vector ^'^^A^. 

Similarly, one can also obtain the following expression 
for the magnetic field Bi in terms of ^'^•'Ai: 

= (3)gim"^^ (3)^^ ^ X (3)^y ^ (2.89) 

with the rotational operator {D x ^^^A) defined in the 
same way as before. 

At this point, one could take the point of view that 
the independent dynamical variables are in fact ^^^A^ = 
h^J ^^^Aj and E^ , with evolution equations given by (|2.37p 



and (|2.88p . and simply define the magnetic field through 
Eq. (|2.89p . therefore ignoring Eq. (|2.38p (which is now 
just a consequence of the definition of i?'). One would 
also find that the magnetic constraint (|2.3ip is now triv- 
ial. We would then have a system of evolution equations 
that is first order in time and second order in space, with 
only one constraint, namely the electric constraint (j2.30|) . 

Of course, one would still have to choose a gauge con- 
dition in order to evolve the scalar potential <i>. One 
possibility would be to take the Lorentz gauge, which is 
given in terms of the 4-vector potential as: 

VaA" = . (2.90) 

This gauge condition can be easily seen to take the fol- 
lowing form in 3+1 language 

9*$ + Cn<P = -Dm (iV (^U™) + NK^ , (2.91) 

with K the trace of the extrinsic curvature, which clearly 
provides us with an evolution equation for $. 

Taking this point of view, however, has one serious 
drawback. One can show that the system of evolution 
equations given by ((07)l . and ([^^ is in fact 

only weakly hyperbolic even in flat space {i.e. it has real 
eigenvalues, but does not have a complete set of eigenvec- 
tors), so that the system is not well-posed. There are cer- 
tainly ways around this, involving defining new auxiliary 
variables and crucially commuting the second covariant 
derivatives of ^^'>A^ that would appear in Eq. (|2.37p when 
we write the magnetic field as (j2.89p (which considerably 
complicates the equation by bringing in a contribution 
from the Riemann tensor), but we will not go into such 
details here (see e.g. Ref. [3a|)- It is enough to say that, 
even though strongly hyperbolic versions of the evolution 
system for '^^^A' and E'' do exist, it is simpler and much 
cleaner to work with the electric and magnetic fields di- 
rectly, and consider the potentials just as auxiliary vari- 
ables when (and if) they are needed. 



III. INITIAL DATA FOR MULTIPLE CHARGED 
BLACK HOLES 

In this Section we will consider the problem of finding 
suitable initial data for multiple charged black holes. For 
simplicity, will concentrate on the case of time-symmetric 
initial data for which the extrinsic curvature vanishes 
Kab = 0. Furthermore, we will also assume that the ini- 
tial magnetic field is zero, in which case the momentum 
constraints are identically satisfied. 

The problem of solving the Einstein-Maxwell con- 
straint equations was studied previously by Bowen 
in [3], where he used a method of images (notably for 
the electric field) to construct a manifold that represents 
two isometric asymptotically flat universes with n throats 
connecting them. In that work Bowen describes a so- 
lution for the electric field that is inversion symmetric 
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through the throats (but which do not arise from a po- 
tential), which then allows one to find numerical initial 
data for the 3-metric that represents n charged black 
holes. Here, however, we will use a different approach 
and will look instead for solutions that are not inversion 
symmetric, but that rather represent a series of throats 
that connect our universe to n distinct asymptotically flat 
universes, more in the spirit of the Brill-Lindquist 37] ini- 
tial data for time-symmetric black holes, or the Brandt- 
Bruegmann p6l | puncture data for spinning or moving 
black holes. 



A. Conformal transformation of the metric and 
electric field 

For time symmetric initial data and vanishing mag- 
netic field, the problem of finding initial data reduces to 
finding a solution of the electric constraint (|2.30[) 

Z?a^" = , (3.1) 

together with the Hamiltonian constraint (|2.79p 

= 167r£ , (3.2) 

with the energy density of the EM given by: 



£ = -5- EaE- 

Stt 



(3.3) 



Let us now assume that the spatial metric hab is con- 
formally fiat, so that we can rewrite it as: 



hab = V' hab 



(3.4) 



with -0 the conformal factor and hab a fiat background 
metric in arbitrary coordinates. The Hamiltonian con- 
straint then reduces to the following elliptic equation for 
the conformal factor 



1 



D^^ + ^ Tp^EaE" ^ , 



(3.5) 



with the Laplacian operator compatible with the 
background metric. 

With the above conformal transformation in mind, let 
us now consider the electric constraint (|3.ip . Notice first 
that Da is the derivative operator associated with the 
physical metric hab- Notice also that, quite generally, for 
any vector we have 



Dav" = Dav" + 6v'' da In 



(3.6) 



with Tp the conformal factor introduced above and D the 
derivative operator associated with the conformal metric 
hab- This implies in particular that 

Da ii^'^v") = i/-" lOav'' + (6 + n)v''da In ip] - (3.7) 



Using this result it is then natural to define the con- 
formally rescaled electric field as 

E'' := ^'^E'' , Ea := ^^Ea , (3.8) 

and the conformally rescaled charge density as 

p ■-= ip^'p . (3.9) 

The electric constraint then reduces to 

Da^" = 47rp , (3.10) 

where now the divergence is calculated with respect to 
the conformal metric. 

In terms of the conformal electric field just defined, the 
Hamiltonian constraint takes the final form: 



D^iP + -i^ EaE'' = 



(3.11) 



In order to find initial data, one must then first solve 
the conformal electric constraint (|3.10[) . and then plug 
in the solution for i?" into the Hamiltonian constraint in 
order to solve for the conformal factor ^p. 

In fact, Eq. (|3.1ip can also be written as 



1 



1 



^D'^P- -{D''^){Dai^) + -EaE^ ={) . (3.12) 

where ip := ip"^ . This equation was considered by 
Bowen [23| for solving the initial data for the single 
charged black hole case (see the analysis below). 



B. Exact initial data multiple charged black holes 
with the same charge-to-mass ratio 

In order to find initial data for multiple charged black 
holes we will first assume that the background metric is 
fiat. Let us introduce a conformal electric potential 
such that 



Ea = -da^p 



(3.13) 



Do notice that the conformal potential p does not coin- 
cide with the physical potential <I> discussed in Sec. Ill Dl 
above, since even in the absence of a vector potential 
equation (|2.88p clearly shows that the relation between 
<& and the physical electric field involves the lapse func- 
tion. 

Using Eq. (|3.13p . the electric constraint can be rewrit- 
ten as 



D'^Lp = -A-Kp 



(3.14) 



From now on we will also assume that we are in a region 
away from any charges, so that p = 0. 

Before attempting to find initial data for multiple 
charged black holes, let us recall for a moment the 
Reissner-Nordstrom analytic static solution for a single 
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charged black hole with mass M and charge Q [ij, llj| , 
for which the conformal electric potential Lp and the con- 
formal factor ip are given by: 



Q 



2r ) 4r2 



1/2 



(3.15) 



where we have assumed that the black hole is centered 
on the origin of the coordinate system r = 0. 

The conformal and physical electric fields for this so- 
lution are purely radial and are given by: 



6' = ^ 



(3.16) 



The fact that the conformal factor -0 above is an exact 
solution of Eq. p. lip for this electric field can be verified 
by direct substitution. 

Since in this case the spacetime is static, Eq. ()2.82p 
provides a linear elliptic equation for the lapse: 



8ttN£ = NE^Ea 



which in terms of the conformal variables reads 

2 - - N - - 

D^N + -{D'N){D,ij) = — EaE'' . 



(3.17) 



(3.18) 



One can now also confirm by direct substitution that the 
lapse given by Eq. (|3.19p below solves Eq. (|3.18p when 
using in turn Eqs. (|3.15p and (|3.16p : 



N : 



(1 + M/2r)(l - M/2r) + Q^/Ar^ 
(l + M/2r)2 -Q2/4r2 



(3.19) 



Notice that Eqs. ((XT5)) . ((XTC)) and ((XTg|l correspond 
to the Reissner-Nordstrom solution in isotropic {i.e. 
conformally flat) coordinates, and not in the standard 
Schwarzschild-like coordinates one flnds in most text 
books It is clear that by taking Q = 0, the above so- 
lution reduces to the Schwarzschild solution in isotropic 
coordinates. 

Based on the form of the conformal factor for the 
Reissner-Nordstrom solution given above, we will now 
propose the following ansatz for the conformal factor in 
the presence of a generic electric potential ip that is a 
solution of the electric constraint (|3.14p : 



(3.20) 



Substituting this back into equation p.l2p we find, after 
some algebra, the following elliptic equation for rj: 



-^^^^^^^-(^-±^) = 0, (3.21) 



where we have already used the fact that D'^(p = 0. 

We can now easily notice a remarkable fact: If we take 
the function rj to be proportional to the electric poten- 
tial ip, that is r] — kip for some constant fc, then equa- 
tion (|3.12p is identically satisfied since in such a case we 
clearly have dm {"n/v) — and iP'rj = (remember that 
ip solves the electric constraint away from charges). 

We will now use this fact to find an exact solution of 
the Hamiltonian constraint for multiple black holes. Let 
us assume that in the conformal space we have a series 
of point charges with values Qi located at the points r^. 
The solution for the potential tp is then clearly 



E 

i=l 



(3.22) 



Let us now choose r/ proportional to ip in the following 
way: 



= kip ~ k 



r^|r-r,| frt2|r-f,| 



(3.23) 



We can now construct a conformal factor that satisfies 
the Hamiltonian constraint as in (|3.20p : 



1+E 




(3.24) 



with Mi = 2kQi, and k an arbitrary constant. This 
solution represents a series of n charged black holes, all of 
which have the same charge-to-mass ratio Qi/Mi — l/2k. 

One could ask at this point how can we know that 
the above solution in fact does represent a series of black 
holes. There are several ways to see that this should be 
so. First, notice that if we take n = 1 this is just the stan- 
dard Reissner-Nordstrom solution for a single charged 
black hole. Second, in the case when all the charges 
vanish the conformal factor p.24p above reduces to the 
well known Brill-Lindquist conformal factor for a series of 
non-charged black holes [37]. Also, if the different points 
fi are very far apart, then close to each of them the con- 
formal factor again reduces essentially to the Reissner- 
Nordstrom solution, so we would indeed have a series of 
n charged black holes. Finally, notice that because of 
the singularities in the conformal factor, as we approach 
each point the areas of spheres centered around that 
point first become smaller and then increase again. That 
is, our initial data is a topological construction with a 
series of wormholes connecting to other asymptotically 
fiat regions. Since the initial data is time symmetric the 
throats of these wormholes (i.e. the minimal surfaces) 
in fact correspond to apparent horizons. Of course, if 
some of the points fi are very close to each other one 
could find common apparent horizons around them, so 
that we actually have fewer black holes with complicated 
internal topologies. As a final comment, notice also that 
due to the tidal forces between the different black holes, 
the throats of the wormholes can not be expected to be 
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spherical, and their precise shape and location should be 
found numerically. 

That an exact solution of the Hamiltonian constraint 
for multiple charged black holes exists at all is a surpris- 
ing result. Notice, however, that we have only found a 
solution for the initial data. In general, one would expect 
this initial configuration to evolve as each black hole re- 
acts to the gravitational and electric fields of the other 
black holes, so that a non-zero extrinsic curvature and 
magnetic field would quickly develop. 

There is in fact one notable exception to this. In order 
to find it we will first rewrite the conformal factor (|3.24p 
in the following way 



V'^ = l + 2fc^- 



Qi 



charge-to-mass ratio. When considering different charge- 
to-mass ratios we have not been able to find a closed form 
solution. On the other hand, finding numerical solutions 
can be done easily enough. In order to do this we will 
first modify our ansatz (|3.20p above in the following way: 



ip = ip"^ — {u + rj) — 



V5 



with 



N 



V = J2 



M,. 



-{2\T-n\ 

N 



(3.28) 

(3.29) 
(3.30) 



(3.25) 



where we have already used the fact that Mi — 2kQi. If 
we now take k — 1/2, which implies Mi — Qi, then the 
conformal factor reduces to: 



^2 



n 

E 

1=1 



Q^ 



(3.26) 



which now corresponds to initial data for a series of ex- 
tremal black holes. Amazingly, one can show that this 
multi-extremal black hole solution turns out to be static, 
that is, the gravitational attraction of all the black holes 
is exactly canceled out by their electrostatic repulsion, so 
that the black holes never move from their initial posi- 
tions. The lapse function for such a static solution is 
given by: 



N = Ijil? 



1+E' 



(3.27) 



which solves Eq. I|3.18p exactly. 

This multi-extremal static solution was first obtained 
by Papapetrou TB] and Majumdar [l^ (see [l^ for a 
review), and was further analyzed by several authors (TtI . 
[ii,[3i,[3i,[43| [el]. This solution is so well known that it 
can even be found in some text books (see e.g. the recent 
book by Carroll 41]). On the other hand, as far as we are 
aware the exact solution of the Hamiltonian constraint 
for non-extremal black holes with equal charge-to-mass 
ratios presented above was not previously known. 



and where m is a function that goes to 1 at infinity, and 
in fact is identically equal to 1 everywhere in the case 
when all charge-to-mass ratios are equal. 

Substituting now Eq. (|3.28p into the Hamiltonian con- 
straint we find the following elliptic equation for u: 



4 

r/ + u — 1 



{t) + u) ( (r/ + uX-^\ D^u 



0, (3.31) 



The above equation needs to be solved numerically for u 
in the case when the charge-to-mass ratios of the different 
black holes are not all equal. 

Before presenting some examples of numerical solu- 
tions for the case of two black holes, it is important to 
investigate the expected behaviour of the function u close 
to each of the "punctures" , that is, close to each of the 
points r — fi. In order to do this, let us now use a sys- 
tem of spherical coordinates (r, 0, </>) adapted to one of 
the black holes. Without loss of generality we will choose 
that black hole as the one identified with the label 1, so 
that ri = 0. Let us now examine the behavior of the 
different terms in equation (|3.3ip for small r. Consider 
first the coefficient of the Laplacian operator: 



Ti := {ri + u)[ {ri + uf - 



(3.32) 



Let us assume for the moment that u is finite at each 
of the punctures. From the form of the functions 77 and 
(p, it is then clear that for small r this term behaves in 
general as 



C. Numerical initial data for multiple charged 
black holes with different charge-to-mass ratios 

The exact multiple black hole solution of the Hamil- 
tonian constraint found in the previous section is only 
valid in the case when all the black holes have the same 



Ti ~ l/r^ 



(3.33) 



Consider now the term with first order derivatives in 
Eq. (IH3TD : 



(3.34) 
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In order to analyze the behaviour of this term for small 
r, we will first expand the derivatives to obtain 



D. Numerical examples for the case of two charged 
black holes 



T2 = 9r„(7? + u) - (77 + U - 1) dm^] 

[(pa'"(7? + u)-(77 + u+l)9"(^] . (3.35) 
We will now rewrite the functions 77 and above as 

V=^+H{r,9,^) , ^^9l + F{r,9,(P) , (3.36) 
2r r 



with H and F given by 

Mr 



2\r-n 



(3.37) 



which are clearly regular functions at r = 0. Substituting 
into T2 and expanding we find, after some algebra, that 
for small r this term behaves as: 



T2 ^ l/r^ 



(3.38) 



Notice that naively one could expect T2 to diverge as 
1/r^, due to the presence of terms of the form ((/? drrj)^, 
but in fact all such terms cancel out and we are left with 
a dominant divergence of order l/r*. 

From the behaviour of Ti and T2 for small r, we then 
find that in order for equation (|3.31[) to be consistent the 
Laplacian of u must behave for small r as: 



1/r . 



(3.39) 



Now, in spherical coordinates the (flat) Laplacian is 
given by 



2 „ 1 



D^u ^d^u+^ dru +-^L^u, (3.40) 



with the angular operator 



L^u := de (sin 9 deu) + -A^- dlu . (3.41) 
smt' o ¥^ 



1 

s\^9 



This implies that in order to have D^u behaving as ex- 
pected for small r we must ask for the function u to have 
a Taylor expansion near the origin of the form: 

u = a + h{9,<t>)r , (3.42) 

with a a constant and h{9, (j)) some regular function of the 
angular coordinates. Notice that b{9, (p) must be non- 
zero for equation (|3.3ip to be consistent, so the above 
expansion implies that u is not regular at the origin (re- 
member that r is a radial coordinate). The function u 
then turns out to be only C° at the punctures, that is, it 
is finite and continuous, but its derivatives are no longer 
continuous. In other words, the function u is in general 
expected to have a kink (i.e. a conical singularity) at 
each of the punctures (we will see in the numerical ex- 
amples below that this is indeed the case). 



We have constructed a simple numerical code to solve 
equation p.3ip for the case of two charged black holes 
with different charge-to-mass ratios. In such a case one 
can locate both black holes along the z axis. The sit- 
uation is then clearly axisymmetric, so the problem is 
effectively two-dimensional. 

As a boundary condition we ask for the function u to 
behave as w = 1 -|- c/r for large r, where r = -^/p^ + 
and c is some constant. In order to eliminate the un- 
known constant this boundary condition is differentiated 
and applied in the following way: 



drU 



1 



(3.43) 



Our code uses cylindrical coordinates (p, z, (j)) instead 
of spherical coordinates {r,9,(j)), so that in practice we 
assume that far away the dependence on the azimuthal 
angle 9 can be ignored, so that one can write 

dpU = (^) drU , d^u = (^) drU . (3.44) 

The final boundary condition on the p boundaries is then 



r\ l-u 

-]dpu= 

P T 



(3.45) 



with an analogous condition on the z boundaries. 

Since here we are mainly interested in showing that the 
solutions for u can be easily found and that they behave 
as expected, we have decided to write a very simple code 
that instead of solving the elliptic equation directly solves 
an associated hyperbolic problem of the form: 



u-1 



1 



(3.46) 



with t a fictitious time parameter. The boundary condi- 
tions are then also modified to outgoing- wave boundary 
conditions of the form: 



dtu - 



P T 



(3.47) 



and analogously for the z boundaries. We choose as ini- 
tial condition u = 1, and evolve the above hyperbolic 
equation until we reach a stationary state with some 
predetermined tolerance. The idea is that the wave-like 
equation above will propagate the residual away through 
the boundaries, and will drive the system to a static so- 
lution which corresponds to the solution of the original 
elliptic problem. 

Spatial derivatives are approximated with standard 
centered second order differences, while for the time in- 
tegration we use a three-step iterative Crank-Nicholson 
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algorithm [42, |43|- The resulting code is rather slow, 
particularly for high resolutions, but has the advantage 
of being both very simple to write (and debug), and also 
extremely robust. A more sophisticated elliptic solver 
should of course be used in order to find highly accurate 
solutions in a short computational time (we certainly do 
not recommend to use our quick-and-dirty "wave-like" 
algorithm for any kind of production runs). 

In order to deal with the problem of divisions by p 
our computational grid staggers the axis and introduces 
a ghost grid point at position p — —Ap/2. To determine 
the value of u at this ghost point we then simply impose 
the parity condition u{—Ap/2) = u{Ap/2). 

Notice that we don't do anything special close to the 
punctures and just take standard centered differences ev- 
erywhere. This clearly affects the order of convergence 
close to the punctures (see numerical examples below). 




Example I: Equal masses and opposite charges 



FIG. 1: Numerical solution for the function u along the z axis 
for the case of equal masses and equal but opposite charges. 



As a first example we will consider the case of 
black holes with equal masses Mi — M2 — 1, and ec 
but opposite charges Qi = —Q2 = 1/2. The puncti 
are located along the z axis at positions zi — — Z2 = 
and the boundaries extend to p = 6, z = ±6. 

As initial guess we choose u = 1 everywhere, and 
then evolve equation p.46p until the magnitude of 
right hand side is everywhere smaller than a fixed to 
ance of e = 10~^. This tolerance is chosen such tha 
is always smaller than the truncation error at the resi 
tions considered here. 

Figure [1] shows the numerical solution for the fund 
u along the z axis, for a resolution of Ap — Az = 0.01 
Notice how the function u behaves as expected on b 
punctures, with very evident kinks. One can also see t 
u has equatorial symmetry even though the charges h 
opposite signs. Of course, had we chosen both chai 
with the same sign the exact solution would have b 
u = 1 everywhere since in that case the charge-to-niM,^,u 
ratios would have been equal. It is also interesting to 
note that the maximum deviation of the function u from 
unity is precisely at the punctures, and this maximum 
deviation is of only about 7%. Figure [2] shows the a 
height map of the same solution in the {p, z) plane. 

Next we show a plot of the convergence in the Hamilto- 
nian constraint along the z axis. Since we are in fact solv- 
ing numerically precisely the Hamiltonian constraint, one 
would expect it to be satisfied to the level of the tolerance 
in the elliptic solver. This is of course true, but in order to 
study convergence we are in fact using a different expres- 
sion for the Hamiltonian constraint. We first reconstruct 
the conformal factor and later evaluate numerically to 
second order the Hamiltonian constraint p. lip written 
as: 




^V+T^5™¥'5"V = 
4^' 



(3.48) 



FIG. 2: Height map in the (p, z) plane of the numerical solu- 
tion for the function u. 



This last expression should not be expected to hold to 
the level of the tolerance in the elliptic solver, but rather 
to the level of numerical truncation error, which is much 
higher (we have chosen a very small tolerance parameter 
in the elliptic solver precisely for this reason). 

Figure [3] shows the logarithm of the absolute value 
of the Hamiltonian constraint evaluated using (|3.48p 
for the five different resolutions Ap — Az ~ 
0.1,0.05,0.025,0.0125,0.00625, with each plot rescaled 
by the corresponding factor expected for second order 
convergence: 1,4, 16, 64, 256 (we in fact show only the re- 
gion close to one of the punctures as the situation is sym- 
metric on the other puncture) . Notice how away from the 
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log(ham)[0.1' 
log(4*ham)[0.05' 
I log(16*ham)[0.025' 
log(64*ham)[0.0125' 
|log(256*ham)[0.00625' 




0.5 1 



1.5 



2.5 



3.5 




FIG. 3: Convergence of the Hamiltonian constraint. We plot 
the logarithm of the absolute value of the Hamiltonian con- 
straint evaluated at 5 different resolutions, with each reso- 
lution rescaled by the corresponding factor in order to show 
second order convergence. 

puncture all plots lie on top of each other, showing nice 
second order convergence. Closer to the puncture, how- 
ever, the Hamiltonian constraint in fact increases with 
higher resolution, but this loss of convergence is lim- 
ited to the 2 or 3 grid points closest to the puncture, 
so that the non-converging region keeps getting smaller 
and smaller with higher resolution. This should not be 
surprising since in the conformal factor ip we have terms 
that diverge as l/r close to each puncture. 



2. Example II: Equal masses and one charge equal to zero 

As a second example we will consider the case of two 
black holes with equal masses Mi — M2 — 1, and one 
charge set equal to zero Qi = 1/2, Q2 = 0. As before, 
the punctures are located along the z axis at positions 
zi = —Z2 = 2 and we use a resolution of Ap — Az — 
0.0125. The results are plotted in Figure H Notice that 
even though Q2 = 0, the function u still has a small kink 
at z — —2. 



3. Example HI: Some more generic cases 

As a more generic case we set up two black holes 
with different masses. Mi = 1,M2 — 0.5, located as 
before at the points zi = —Z2 = 2. The first black 
hole has a charge of Qi — 1/2, while for the charge 
of the second black hole we consider a series of values: 
Q2 = -0.25,-0.1,0.0,0.1,0.25. The results are plotted 
in Figure O 

What seems to happen as Qi > is kept fixed and 



FIG. 4: Similar to Fig.[T] but for the case of equal masses and 
one charge equal to zero. 
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Q2=0.00 






Q2=+0.10 
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Q2=+0.25 - - 
1 1 1 





-6 -4 -2 2 4 6 

Z 

FIG. 5: Similar to Fig. [T]for a more generic case in which the 
back hole at z = 2 has Mi — 1, Qi = 1/2, while for the black 
hole at 2: = —2 we take M2 = 1/2 and a series of values for 
the charge: Q2 = -0.25, -0.1, 0, 0.1, 0.25. 

Q2 changes is that, for Q2 < 0, the kinks in u are point 
down at both punctures. Then, as Q2 approaches zero 
the kink at that puncture becomes much smaller but still 
keeps pointing down. If Q2 goes through zero and be- 
comes positive, then the kink at Z2 first disappears and 
later starts growing again in the opposite direction. As 
Q2 keeps increasing towards a value where both charge- 
to-mass ratios are equal both kinks become smaller with 
opposite directions, until the function u becomes 1 ev- 
erywhere. 

By trial and error we have in fact found a situation for 
which the kink at z = —2 seems to essentially disappear. 
Figure [S] shows the numerical solution for u in the case 
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0.995 - 



0.99 - 



0.985 - 

0.98 I 1 1 1 1 1 

-6 -4 -2 2 4 

z 

FIG. 6: Solution for u in the case when Afi = 1,M2 — 
0.5, Qi = 0.5, Q2 = 0.01475. Notice that the kink at z = -2 
is essentially gone. 



when Ml = 1, M2 = 0.5, Qi = 0.5, Q2 = 0.01475. 

The special value of Q2 for which the kink disappears 
can in fact be estimated analytically. When we examined 
the behaviour of the function u near the punctures in 
Section IIII CI above, we mentioned the fact the source 
term of equation p.3ip . which we called T2, behaves as 
1/r^ near each puncture. Being somewhat more precise 
one can show that for the case of two black holes this 
term behaves close to the puncture at r2 as 

= 4^ [— d '^^^"j 

X I — 4(32 - 2(92Au 1 

+ 0(l/r3) , (3.49) 

with r — |r — r2|, and where d is the distance between the 
punctures and Au = 14(^2) — 1- Notice that there are two 
ways in which the term T2 can vanish, corresponding to a 
solution u that is regular at the puncture. One possibility 
is to have u — 1, so that Am = 0, and equal charge-to- 
mass ratios so that Q1M2 = Q2M1, in which case the first 
term above vanishes. This is clearly the exact solution we 
already described. But a second possibility is to assume 
that u is very close to 1 at the puncture so that Au ^ 1, 
and to ask for Q2 ^ Q1M2I {Mi + 4(i), in which case the 
second term above will almost vanish (T2 will in fact not 
vanish exactly for this value of Q2 since at the puncture 
Au is not exactly zero, but one can expect T2 to vanish 
for a value of Q2 that is very close to this one). For 
the values in our example we have Mi = 1, M2 = 0.5, 
Qi = 0.5 and d = 4, so that Q2 - 1/68 - 0.0147, which 
is remarkably close to the empirical value found above. 



We have considered the Einstein-Maxwell system hav- 

- ing two goals in mind. The first goal consisted in recast- 
ing the covariant Maxwell equations in a curved space- 

_ time as an initial value problem along the lines of the 
usual 3+1 formalism of general relativity, choosing the 
magnetic and electric fields as independent variables. 

- This led to a set of two constraint equations for the elec- 
tric and magnetic fields, plus a set of two evolution equa- 
tions for those fields. The evolution equations are hy- 

- perbolic with the propagation speeds depending on the 
lapse and shift. 

The second goal was to construct initial data satisfying 
6 the gravitational and electromagnetic constraints which 
represent momentarily static charged black holes. In or- 
der to achieve this goal we assumed a moment of time 
symmetry with vanishing extrinsic curvature and mag- 
netic field. For the case of two black holes, this initial 
data will serve to analyze a head-on collision. We found 
that for black holes having an equal charge-to-mass ratio 
it is possible to find an analytic solution for the con- 
straints. However, when this condition is dropped, we 
present instead numerical solutions. We studied the be- 
havior of such solutions with several values of the free 
parameters {i.e. different masses and charges), showing 
extremal cases. 

A much more realistic initial data will consist in relax- 
ing the moment of time symmetry condition {i.e. aban- 
don the initial condition Kij = 0), but keeping a null 
initial magnetic field and analyze the analogous of the 
Bowen-York initial data. In this case, one will need to 
solve a much more complex Hamiltonian constraint for 
the conformal factor. An elliptic solver will be required 
a fortiori to solve this constraint. This is an issue that 
is worth analyzing in the future. 

Using the initial data presented here we plan in a near 
future to study the evolution of the Einstein-Maxwell sys- 
tem generated by the collision of two charged black holes, 
and analyze the emission of both gravitational and elec- 
tromagnetic radiation. We believe that the interplay be- 
tween gravity and electromagnetism will give rise to very 
interesting and unusual dynamics due to the possibility 
of repulsion between charges. This is an issue that has 
not been scrutinized thus far in numerical relativity. 
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APPENDIX 

For the sake of clarity, we will give here some details 
about our derivation of the 3+1 Maxwell equations. We 
start with Eq. (|2.ip . The projection of this equation onto 
ria reads: 



(1) 



where one must remember p :— —ribj^. At this point 
there are several possible ways to proceed. Let us first 
write 



(2) 



where in the last step we used Eq. (|2.18p . 

Now, concerning the term Vai?°, one way to relate it 
with 3+1 quantities is by using the following identity for 
any 4-vector field 



da {^gV'') 
da (NVhV 



(3) 



where we have used the fact that g 
have 



-N^h. We then 



V^V^ = V'^dah-iN +^da [VhV^ 



(4) 



For the particular case where = E"-, it turns out 
that E'^dalnN = E'^/i^'^ValniV = E^DalnN, where the 
fact that E'^ = E^h^^ is a consequence of the fact that 
is a 3-vector. Moreover, since E* = the term 
(l/Vh) da{VhE°-) reduces to the 3-divergence DaE°-. 
We then have 



VaE" = DaE"- + Eaa" , 



(5) 



where we used the following identity for the acceleration 
of the Eulerian observer Of, :— n°'V arib = Di,\nN. We 
point out that one could have obtained the same result 
from the definition of the 3-covariant divergence: 

DaE" = h^^h^VcE'' ^ h.^VcE'' = VbE'' + n^nbVcE'' 



= VbE" 



n^Wc (ubE'') - n'-E^Vcnb 



= VbF^-E^ab 



(6) 



where we used the fact that UbE^ = since is by def- 
inition orthogonal to n". From this equation one obtains 
the same identity Q . Using this result Eq. ^ now reads 



Ub^aF"'' = DaE'' + Eaa" - F^'^am 



(7) 



We now need to prove that, in fact, F'^^V anb — Eaa"' ■ 
For this we use Eq. (|2.20p to obtain 

F'^^'VaUb = (^'^^"'Vanf, + n'^E^'Vanb - E'^n^'Vanb. (8) 



The second term in the last equation is precisely 
■nf^E^V anb — E^'ttb, while the third vanishes identically 
since n^nb = — 1, and therefore Va (n^Ub) — 2n^Wanb = 
0. Finally, we show that the first term vanishes by using 
the definition of the extrinsic curvature, Eq. (|2.8p : 



^""^F-'Kab = - '^^^F'^'h^h.^Vcfia = - '^^^F-^VcUd. (9) 

This term vanishes identically because '^^^F"''' is antisym- 
metric while Kab is symmetric. We have thus proved that 
F°''^Vanb = Eaa°-. Using this result in Eq. ([7]) allows us 
to conclude that 



(10) 



Following an analogous procedure, one can now also 
obtain Eq (|2.3ip by projecting Eq. (|2.2p , except that one 
now uses Eq. (|2.29p instead of (|2.20p . In this case, there 
is no magnetic charge. 



We proceed now to obtain Eq. (|2.37p by projecting 
2?T|) onto St using h°■^^■. 

h'^^VcF"'' ^ ^An^^^f . (11) 

Using Eq. (|2.20p and expanding we obtain 

/i°bVc ''^^F"'' - E'=h''t,Vcn'' + E'^Vcn'' + h'^^n'^VcE'' 



h^.Vc '■^^F"'' + E''K.'' - E^K + h^.n^^/cE^ 



~47r j 



(3) Aa 



(12) 



where we used E^h^^X/^n'' = E'^h/h''^\/ cu'' = -E'^Kf-, 
and K = -\/ [cf. Eqs. (gH) and ([fill)) ]. 

Now, from the definition of the 3-covariant derivative 
applied to a 3-tensor field we have [cf. Eq. (|2.1ip ] 

- h'^fV/^'^F'^f- ^^^F^'^Vdinen") 
In this way Eq. becomes 



h''j,n''VcE^ + De (3)^ca _^ (3)^ 



ba 



ab 



(3) ,a 



(14) 



The first term above can be rewritten using the following 
identities 



h^uTi'^WcE'' + K'^.E" 



(15) 



We then conclude that 



h\CnE^ + Df^F"'' ^^^^F^^ab 



'E^K = -At: '■•^'j 



(3),-a 



(16) 
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Finally, using the fact that (3)^ca jg antisymmetric, 
together with Eqs. (|2.27p and (|2.25p . one can write 

= ^ SbB, = - . (17) 

Equation (fTB|) then reads 

-^E^'K ^ -iTT^^T , (18) 

where we used again Eq. (|2.27p to rewrite the third term. 
Taking the notation given by Eqs. ()2.34|) and (I2.35p . one 
finally recovers Eq. (|2.33p . 

In order to find Eq. (|2.37p from Eq. (fT8|) . one has to 
remember that the Lie derivative can be written in terms 
of ordinary derivatives as follows 

h^^CnE^ = h''^{n''dcE''~E''dcn'') . (19) 

Now, since we are only interesting in the spatial com- 
ponents of Eq. (fTB]) . using (|2.5p together with = 
(-7V, 0, 0, 0), n° = (1/7V, N'/N) one can easily show that 

h\C^E^ = ^ d^E^ + ^ d,E^ ~ ^ d,N^ . (20) 

Substituting this last result into the spatial components 
of HHD leads directly to Eq. (|07)) . 



In a similar fashion one can obtain Eq. (|2.38p by pro- 
jecting (|2.2p onto S(. In fact, this just amounts to using 
the duality relations Ea Ba and Ba —Ea in (|2.37p 
to obtain Eq. (gSSl)- 

Finally, we proceed to derive Eq. (|2.40p . This can be 
easily done by using the orthogonal decomposition of the 
electric-current 4- vector [cf. Eq. (|2.6p ]: 

f = ^^^f+pn" , (21) 

being p :— —Uaj"" the charge density measured by the 
Eulerian observer, and := h'^^j'^. This decomposi- 

tion, when inserted into Eq. p.39p and making use of 
(PrU)) . leads directly to 

n'^VaP -pK + Va ^^^f = . (22) 

Using now Eq. dSl) for (3) instead of E°- (the expression 
is valid for any 3- vector), one obtains 

CnP + Da ^^\f + ^^^faa - pK = , (23) 

where 

£„p EE n'^WaP = l/N {dtp + Wd,p) . (24) 
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